##########################################################################################

library(data.table)
library(optparse)
library(Seurat)
library(ArchR)
library(ggthemes)
library(dplyr)
library(pheatmap)
library(ggplotify)
library(cowplot)

##########################################################################################
option_list <- list(
    make_option(c("--rna_file"), type = "character"),
    make_option(c("--comine_data_all_file"), type = "character"),
    make_option(c("--known_gene_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){

    ## 单细胞表达文件
    rna_file <- "~/20231121_singleMuti/results/qc_atac/testis_combined.annotationCellType.qc.Rdata"

    ## 所有的细胞的,计算maxdelt
    comine_data_all_file <- "~/20231121_singleMuti/results/qc_atac/testis_combined_peak.combineRNA.qc.Rdata"

    ## 认为关键的TF
    known_gene_file <- "~/20231121_singleMuti/config/SCOS_known_genes.csv"

    ## 输出
    out_path <- paste0("~/20231121_singleMuti/results/pheatmap_knowngene")

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

rna_file <- opt$rna_file
comine_data_all_file <- opt$comine_data_all_file
known_gene_file <- opt$known_gene_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################
## 导入数据
a <- load(rna_file)
## scrnat
DefaultAssay(scrnat) <- "RNA"

b <- load(comine_data_all_file)
## testis_combined_peak_combineRNA

## 已报道的TF
known_gene <- unique(data.frame(fread(known_gene_file)))

###########################################################################################
## 细胞顺序
cell_order <- c("SSC" , "Differenting&Differented SPG" , "Leptotene" ,
    "Zygotene" , "Patchytene" , "Diplotene" , 
    "Early stage of spermatids" , "Round&ElongateS.tids" , "Sperm" ,
    "Sertoli cells" , "Leydig cells" , "Myoid cells" ,
    "Pericytes" , "Macrophages" , "Endothelial cells" ,
    "NKT cells"
    )

###########################################################################################
## ATAC的开放矩阵
GSM_se <- getMatrixFromProject(testis_combined_peak_combineRNA, useMatrix="GeneScoreMatrix")
GSM_mat <- assays(GSM_se)$GeneScoreMatrix
rownames(GSM_mat) <- rowData(GSM_se)$name

## 转化为数值矩阵
GSM_mat_num <- apply(GSM_mat , 1 , as.numeric)
rownames(GSM_mat_num) <- colnames(GSM_mat)
GSM_mat_num <- t(GSM_mat_num)

## atac的细胞列表
atac_cell_data <- testis_combined_peak_combineRNA@cellColData$cell_type

###########################################################################################
## 构造输入基因列表
all_gene <- known_gene$Human_Symbol
tf_gene <- subset(known_gene,Human_TF.TC!="")$Human_Symbol
non_tf_gene <- subset(known_gene,Human_TF.TC=="")$Human_Symbol

known_gene_list <- list(
    all = all_gene ,
    tf = tf_gene ,
    non_tf = non_tf_gene
    )

for( type in names(known_gene_list) ){

    print(type)
    use_gene <- known_gene_list[[type]]

    ## 构造表达矩阵
    sco_exp <- sapply(unique(scrnat$cell_type),function(x){
        sapply(unique(use_gene),function(y){
            mean(as.numeric(as.vector(scrnat@assays$RNA@data[y,which(scrnat$cell_type==x)])))
        })
    })

    ## 构造atac矩阵
    sco_atac <- sapply(unique(atac_cell_data),function(x){
        sapply(unique(use_gene),function(y){
            mean(as.numeric(as.vector(GSM_mat_num[y,which(atac_cell_data==x)])))
        })
    })

    #聚类热图不可以存在NA值/整行为0
    sco_exp <- sco_exp[which(rowSums(sco_exp)>0),]
    sco_exp <- sco_exp[,cell_order]
    sco_atac <- sco_atac[which(rowSums(sco_atac)>0),]
    sco_atac <- sco_atac[,cell_order]

    ## 基因顺序
    gene_order <- use_gene[use_gene %in% rownames(sco_exp)]

    ############################################
    ## 基因按照聚类结果
    ## 表达矩阵
    p <- pheatmap(t(sco_exp),scale = "column",show_rownames = T ,show_colnames = T, 
             #cutree_cols=7,cutree_rows = 7,
             cluster_rows = F,cluster_cols = T,clustering_method = "ward.D2",
             color = colorRampPalette(c("blue", "white", "firebrick3"))(100),
             cellwidth = 10, cellheight = 10)
    ## 构建ATAC
    p1 <- pheatmap(t(sco_atac),scale = "column",show_rownames = T ,show_colnames = T, 
             #cutree_cols=7,cutree_rows = 7,
             cluster_rows = F,cluster_cols = T,clustering_method = "ward.D2",
             color = colorRampPalette(c("blue", "white", "firebrick3"))(100),
             cellwidth = 10, cellheight = 10)

    ## 按照基因表达聚类
    p <- as.ggplot(p)
    p1 <- as.ggplot(p1)
    out_file <- paste0( out_path , "/pheatmap_" , type , ".cluster.pdf" )
    pdf(out_file , height = 8 , width = 15)
    print(plot_grid(p, p1, labels = c('RNA', 'ATAC'), label_size = 12 , ncol = 1))
    dev.off()


    ############################################
    ## 基因不聚类
    ## 表达矩阵
    sco_exp <- sco_exp[gene_order,]
    p <- pheatmap(t(sco_exp),scale = "column",show_rownames = T ,show_colnames = T, 
        cluster_rows = FALSE, cluster_cols = FALSE, 
        row_order = cell_order, col_order = gene_order ,
        color = colorRampPalette(c("blue", "white", "firebrick3"))(100),
        cellwidth = 10, cellheight = 10)

    ## 构建ATAC
    sco_atac <- sco_atac[gene_order,]
    p1 <- pheatmap(t(sco_atac),scale = "column",show_rownames = T ,show_colnames = T, 
        cluster_rows = FALSE, cluster_cols = FALSE, 
        row_order = cell_order, col_order = gene_order ,
        color = colorRampPalette(c("blue", "white", "firebrick3"))(100),
        cellwidth = 10, cellheight = 10)

    p <- as.ggplot(p)
    p1 <- as.ggplot(p1)
    out_file <- paste0( out_path , "/pheatmap_" , type , ".noncluster.pdf" )
    pdf(out_file , height = 8 , width = 15)
    print(plot_grid(p, p1, labels = c('RNA', 'ATAC'), label_size = 12 , ncol = 1))
    dev.off()
}


